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The  report  contains  the  asymptotic  efficiencies  of  some  candidate 
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measurably  less  effort  than  solving  the  likelihood  equations.  They  include 
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I.  Introduction 


The  need  for  readily  computed  parameter  estimates  is  great. 

Maximum  likelihood  estimators  are  known  to  be  asymptotically  efficient,  but 
often  they  are  very  hard  to  find.  The  most  popular  alternative  is  the  method 
of  moments  which  usually  yields  readily  computed  estimates,  but  sometimes 
these  estimates  are  not  very  efficient.  This  report  looks  at  the  efficiency 
of  the  method  of  moments  and  of  some  similarly  constructed  'quasi-method 
of  moment 'estimators. 

The  basic  idea  is  to  select  a system  of  estimating  equations  which 
equates  various  statistics  to  their  expected  values.  The  method  of  moments 
does  this  for  sample  moments  of  order  one,  two,  etc.  We  propose  to  consider 
also  moments  of  order  zero,  minus  one,  and  perhaps  other  functions.  The 
examination  of  the  consequent  efficiencies  may  aid  in  the  building  of 
intuition  so  that  a wiser  selection  of  estimating  functions  can  be  made  in 
new  situations  when  they  appear. 

Moments  of  order  zero  and  minus  one  require  positive  data.  The  former 

is  the  geometric  mean  and  the  latter  is  the  harmonic  mean.  They  form  part  of 

a general  family,  f(r)  = (-^  of  nondecreasing  means  (where  x^,...,x^ 

fonns  the  data).  When  dealing  with  the  harmonic  mean,  we  will  be  setting 

-S  — (=  [f(-l)]"^)  equal  to  E(l/X)  because,  for  the  parent  populations 
n Xj 

treated  here,  the  latter  value  is  easy  to  obtain.  Similarly  when  dealing 
with  the  geometric  mean,  we  will  be  setting  ^ ^ In  x^  (=  In  f(0))  equal  to 
E(ln  X).  In  this  form,  the  geometric  mean  appears  in  many  maximum  likeli- 
hood systems.  This  suggests  that  alternative  quantities,  that  are  close  to 
means  with  r = 0,  might  be  profitably  exploited.  The  results  give  some 
support  to  this  idea. 

The  structure  of  the  efficiency  computations  utilizes  the  theoretical 

work  presented  in  a companion  report  [12].  The  pertinent  material  is 
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summarized  in  Section  II  of  the  present  report.  Section  III  contains 
applications  of  the  idea  to  two,  one-parameter  parent  populations,  Poisson 
and  symmetric  beta.  Some  speculations  for  interpreting  these  results  are 
made.  Section  IV  is  devoted  to  two  parameter  settings;  gamma,  negative 
binomial,  and  beta. 

Much  computational  work  is  necessary  to  support  this  development  and 
the  details  are  relegated  to  the  appendices.  Appendix  A contains  general 
relationships  among  the  required  population  moments.  Specialization  to 
Poisson,  gamma,  negative  binomial,  and  beta  is  contained  in  Appendices 
B,  C,  D,  E (resp.).  Computations  are  performed  by  APL  programs  and  these 
are  included  too,  as  Appendix  F. 

II.  Methodology 

It  is  assumed  that  the  estimating  equations  have  the  form 


g(x,0)  = 0 


(2.1) 


I 


where  x = (Xj^,...,x^)  is  the  data  vector  of  a random  sample  from  the  specified 
parent  population,  9'  = (9j^,...,9p)  is  the  parameter  point  which  belongs  to 
an  open  subset  of  p-dimensional  space,  and  g'  = , . . . ,gp) . Primes  denote 

transpose.  In  order  for  the  system  (2.1)  to  have  a unique  solution  6,  it 

is  necessary  (by  the  Implicit  Function  Theorem  [16])  that  the  Jacobian 
5gj^f-tg„  I 


0. 


39i,...,9p  I 

The  following  structural  assumptions  are  made: 

(i)  Each  gj(x, 9)  for  j = l,...,p  is  a symmetric  function  of  x,  i.e. 
is  Invariant  under  permutations  of  Xj^,...,x^. 


(ii)  E{gj(X,6)}  = 0 for  j = 

(ill)  Var{g,(X,9)}  = Jc,(0)  + o(-)  for  j = 
j n J n 

(i^)  the  gj(x,e)  have  bounded  continuous  partial  derivatives  with 

respect  to  9^^ 9^  for  j » l,...,p. 

(v)  9,  which  is  the  solution  of  (2.1)  is  consistent  and  asymptotically 

normal. 


Assumption  (i)  is  modest  and  expected  of  any  reasonable  estimating 
scheme.  The  meeting  of  assumptions  (ii)  and  (iii)  is  a matter  of  scaling 
and  arrangement.  Assumption  (iv)  is  needed  so  that  the  asymptotic  covariance 
matrix  is  well  behaved.  Assumption  (v)  is  always  desirable  and  convenient 
since  it  implies  that  the  estimators  are  asymptotically  unbiased  and  the 
ellipsoid  of  concentration  [4]  is  characterized  in  terms  of  the  inverse  of 
the  asymptotic  covariance  matrix.  Efficiency  computations  are  based  on  its 
determinant. 

Equipment  for  verifying  (v)  is  contained  in  [9].  There,  the  functions 

1.  rll 

gj  are  averages  of  the  form  — gj(x^,6)  for  j = l,...,p  and  this 

structure  is  consonant  with  the  present  set  of  assumptions.  All  of  the 

cases  treated  in  this  report  have  this  structure. 

Much  license  is  taken  in  what  follows.  The  purist  is  referred  to  [9]. 

Let  A(x,9)  be  the  p by  p matrix  of  partial  derivatives  {3g./39,  }. 

J K 

Assume  that  the  elements  behave  as  in  (2.2) 

Ajj^(X,9)  ^ E{9gj/99j^}  (2.2) 

as  n The  resulting  limiting  matrix  is  denoted  by  A or  A(9)  . 

The  assumptions  allow  the  first  order  expansion 
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g(x,0)  = g(x,6)  + A(x,0  + p(0-0))(0-0)  (2.3) 

where  p is  a diagonal  matrix  of  random  numbers  belonging  to  the  interval 
[0,1].  Since  the  system  is  soluble,  g(x,0)  = 0 and  we  can  write 

(0-0)  = -A  ^(x,  0 + p(0-0)  g(x,0))  (2.4) 

since  the  continuity  of  A implies  that  of  A ^ and  of  g.  Letting  the 
asymptotic  covariance  matrices,  as  n be  defined  by 

M = limit  nE(0-0)(0-0) ' (2.5) 

C = limit  nE{g(X,0)  g'(X,0)}  (2.6) 

it  follows  from  (2.2)  and  (2.4)  that 

M = a"^  C(A’)"^  (2.7) 


The  method  of  maximum  likelihood  fits  into  this  setting.  The  parent 
population  has  density  f(x,0),  and  (2.1)  takes  the  form 


nr  n 


n 3in  f (x. ,0) 


i=l 


90 


= 0 


for  r = 1, . . . ,p 


(2.8) 


Assumption  (ii)  requires  the  regularity  conditions  [10]  as  does  the  relation- 
ship 

« w 1 

' r k ’ 

where  A is  the  Information  matrix.  Using  (2.2)  and  (2.6)  it  is  seen  that 
both  A and  C are  equal  to  A and  hence 


follows  from  (2.7). 
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(2.10) 


;.,v  ».  •: , j; ; 


Now  suppose  that  only  the  first  q likelihood  equations  (2.8)  are 
used  in  the  system  g = 0.  Let  us  denote  this  subset  by  p = 0,  and  the 
remaining  p-q  equations  by  h = 0.  Thus  in  partitioned  form  (2.1)  becomes 

g = {^  = 0 (2.11) 

h 


Proceeding  formally,  (2.6)  becomes 


( E(pp') 
C = limit  n j 

E(ph')  1 

' j ^11 

^ii\ 

i 

(2.12) 

( E(h'u) 

E(hh')  ] 

1 Si 

^22 ' 

The  information  matrix  can  be  partitioned 

likewise 

1 ''ll 

S2  I 

( 

(2.13) 

1 

1 '^21 

A 22  ) 

Further,  define  a p-q  by 

q matrix 

g^j_  = (E(ohj/90j^)}, 

j = 1 

-.•••.P-q;  k 

= 1 

(2.14) 

and  a 

p-q  order  square  matrix 

g^2  = (EOhjaej^)}, 

j = 1 

■ . • • • »p— q ; k 

= q+i 

(2.15) 

It  is 

shown  in  [12,  Sec.  IV],  that 

! 

I "'ll 

”®21  ) 

C = ' 

f 

^ ~®21 

{ 

^22  ^ 

(2.16) 

where 

is  as  (2.12)  , 

I 

, "'‘11 

~^12  ) 

-'■I 

§21 

1 

®22  ' 

(2.17) 
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^11  ^12 


A2i  H 


(2.18) 


where 


^ ^21^11^12  “ ^21^12®22 


^^21^12^22^  ' ■'■  ^22^2^22 


^11  ~ ■ %2^22®21^ 


^2  ^11^21^22 


(2.19) 


^22  “ ^^22  " ®21^11®n^ 


Because  of  (2.7),  the  determinant  of  (2.18)  is 


Im~^I  = |a|"/|c| 


(2.20) 


and  it  Is  useful  to  draw  attention  to  its  computation.  Using  partitioned 
form,  see  [6,  p.  165],  its  ingredients  have  determinants 


|ci  10^2  - S21^11®12' 


1^11 1 I §22  “ ®21^11^12l 


and  this  is  useful  in  computing  the  asymptotic  efficiency,  [12] 


(2.21) 


Eff(0)  = |m"  ]/|Ai 


(2.22) 


In  the  special  case  of  p = 2 and  q = 1,  (2.19)  reduces  to 


1 


..W 


V • 


(2.23) 


“ “ JcJ  ^^12^2  ■ ^^12^21^22  ^11®22^ 


and  the  determinant  of  M ^ reduces  to  the  especially  convenient  form 


jj^-lj  ^ ^^11®22  " ^12®22^ 


(2.24) 


^11*^22  ®21 


and  the  denomination  of  (2.24)  is  |c| 


III.  Single  Parameter  Settings 
Poisson. 

The  Poisson  random  variable  X has  density 

f(x;X)  = e'^X^/x!  for  x=0,l,..-  (3.1) 

and  the  derivative  of  the  log  likelihood  is 

S = f -1  (3.2) 

It  is  well-known  that 

A = E(S^)  = j (3.3) 

and  the  sample  mean  is  the  minimum  variance  unbiased  estimator  for  all  sample 
sizes.  It  is  also  the  maximum  likelihood  estimate  and  the  method  of  moments 
estimate. 

Let  us  look  further  solely  for  academic  purposes.  Since  the  Poisson 

has  the  property  that  its  mean  is  equal  to  its  variance  it  follows  that  the 
2 

sample  variance  s is  also  a "moment"  estimate  of  \.  Moreover,  the  idea 

,|i 

I 


can  be  extended.  There  are  many  other  statistics  that  can  be  equated  to 
their  expected  values  and  the  resulting  equations  solved  uniquely  for  X. 
One  other  will  be  considered  here,  namely,  the  averages  of  reciprocals  of 


% 

% 


I 


the  {1  + The  one  is  added  as  a convenient  device  for  avoiding 

division  by  zero.  The  result  will  be  called  the  harmonic  mean  estimator. 
2 

Since  s is  directly  an  estimate  of  X,  its  asymptotic  efficiency 
is  quickly  and  easily  expressed.  Using  (2.22),  (3.3),  and  (B.5)  we  have 

Eff(s^)  = ^ (3.4) 

X + 2X 


The  harmonic  mean  estimator  is  characterized  by  the  equation  (see  (B.7)) 


where 


(3.5) 


1 


1 + X. 


(3.6) 


Equation  (3.5)  is  in  the  form  g = 0 (i.e.  (2.1))  and  satisfies  the  assump- 

tions . The  resulting  estimate,  X , can  be  computed  using  the  iteration 
function  that  falls  naturally  out  of  (3.5),  namely, 


X 


r+1 


(3.7) 


and  ^ initializations  Xq  > 0,  but  convergence  can  be 

quite  slow  if  a poor  X^  is  chosen.  Then  asymptotic  variance  C (see 
(2.6))  of  (3.5)  is  the  variance  of  (1  + X)  ^ which  can  be  computed  from 
(B.7)  and  (B.8).  The  quantity  A of  (2.2)  is  obtained  by  differentiating 
(3.5)  with  respect  to  X, 


I 


1 


8 


(3.8) 


Using  (2.22),  (2.7)  we  have 


Eff(X  ) 


A^X/C,  or  irore  explicitly 


Eff(X*)  = ^ (1  - e“^  - Xe"^)^/Var(j^)  (3.9) 

X 

The  efficiences  of  the  two  estimators  are  compared  in  Figure  3.1. 

Of  course  there  is  no  point  in  using  any  estimator  other  than  X = x,  but 
the  graph  suggests  that  the  harmonic  mean  may  be  profitably  used  in  other 
problems  in  which  a choice  must  be  exercised. 


Figure  3.1 

Asymptotic  Efficiency 
( Poisson  ) 


r 1 

Just  for  fun,  let  us  compare  the  values  of  the  three  estimators  x, 

2 * 

s , X when  applied  to  some  famous  data.  First,  the  number  of  deaths  due 
to  mule  kick  in  200  Prussian  army  corps  years.  The  frequency  counts  are 
109,  65,  22,  3,  1 for  variable  values  0 thru  4 [7,  p.  109].  The  estimators 
are 

X = 0.61  , = .6109  , X*  = .6093 

and  agreement  is  rather  good.  Second  (Rutherford  Chadwick  data) , the 
number  of  radioactive  disintegration  in  2608  time  intervals  of  7.5  seconds 
; each  [4,  p.  150].  This  time  we  have 

i X = 3.867  , = 3.633  , X*  = 3.886  . 

The  sample  size  is  much  larger  but  the  value  of  X is  in  a range  of  lower 
2 

efficiency  for  s . Also,  the  radioactive  decay  is  more  properly  modeled 

2 

with  a pure  death  process  and  this  may  help  explain  the  smaller  value  of  s . 

One  final  comment.  The  convergence  of  the  iteration  function  (3.7) 
has  importance  in  finding  the  maximum  likelihood  estimate  for  X from  a 
Poisson  population  in  which  the  zero  values  have  been  truncated.  That  is, 
a trial  that  produces  no  counts  does  not  come  to  the  experimenter's 
attention  [10,  p.  3-13].  The  density  is 

-X  X 

f(x,X)  = TTT  , X = 1,2,  ... 

1 - e ^ 


and  the  maximum  likelihood  equation  is 


Symmetric  Beta 


A Beta  random  variable  with  equal  parameter  values  has  density 


f (x,a) 


X® (l-x)*^  for  Of^x<_l,  0<  a 


(r(o)]‘ 


(3.10) 


Using  the  psi  function  [1], 


d In  r(g) 
da 


(3.11) 


one  can  write  the  derivative  of  the  log  density  as 


= 2i(,(2a)  - 2->(a)  + ln(x(l-x)) 


(3.12) 


A = 2ip'(a)  - 4ij;’(2a) 


(3.13) 


..I.  ■— 

by  (2.9).  Let  In  x = — In  x^  and  express  the  maximum  likelihood 


equation  as 


ip(a)  ~ ip(2a)  = j {In  x + In(l-x) ) 


(3.14) 


so  a,  the  maximum  likelihood  estimate,  is  a function  of  the  geometric 
mean  of  x and  1-x.  Also,  it  is  difficult  to  find. 

Let  us  turn  to  the  method  of  moments.  Because  of  symmetry,  see  (E.ll), 


(E.12) , 


u - 2 ’ 


2 1 


4(2a  + 1) 


Tlius  X cannot  be  used  and  we  turn  to  the  second  moment.  The  form  g = 0 

becomes 


g2  _ 1 ^ Q 

4(2a  +1) 


(3.15) 


ind  the  estimator  can  be  found  explicitly 


1 1 

,2  " 2 


(3.16) 


and  using  (2.2),  (2.7)  produces 


n Var(oi)  a 4(2a  + 1)^  Var(s^) 


(3.17) 


The  right  side  of  (3.17)  requires  that  (E.2)  be  used  with  (A. 2).  The  result 
is  not  compact  and  is  produced  only  by  computer  program. 


A harmonic  mean  option  is  available  since 


a - 1 ^'■1-X^ 


(3.18) 


(see  (E.13))  for  a > 1,  and  the  variances  are  finite  if  a > 2.  Easy 
calculations,  exploiting  (E.14)  and  (E.15),  show  that  both  y and  z, 
where 


T ^ T 

y » i ^ i 

n 1 Xi 


1^1 
n ? 1-x,  ’ 


(3.19) 


should  be  used  instead  of  either  one  alone  in  order  to  reduce  variability 
We  select  the  form  g = 0 to  be 


(a  - l)(y  + z)  - (4a  - 2)  = 0 


(3.20) 


and  solve  explicitly 


* y + z - 2 

I S=  i ■ 

y + z - 4 


(3.21) 


Application  of  (2.2),  (2.7),  (E.14),  and  (E.15)  produces 


r,  var(,*)  . ~ - O' 

a - 2 


for  a > 2 


(3.22) 


. .‘.fc  •_  V 


Figure  3.  2 

Asymptotic  Efficiency 
( Symmetric  Beta  ) 


The  efficiencies  of  these  two  estimates  are  compared  in  Figure  3.2. 

The  second  moment  estimator,  o,  is  uniformly  better  than  the  harmonic  mean 
* 

estimator,  a , and  this  result  needs  interpretation  in  the  light  of  the 

success  in  the  Poisson  case.  This  time  the  averages  of  X ^ and  (1-X)  ^ 

were  used  instead  of  (1+X)  the  latter  being  difficult  to  manage  with 

the  beta  distribution,  and  the  parameter  must  be  at  least  two  in  order 

for  the  variance  to  exist.  Also  the  positive  sample  space  is  0 < x < 1, 

which  entails  large  and  variable  values  for  the  {x^^},  and  this  may  explain 

* 

the  lack  of  success.  The  estimator  ot  may  still  have  some  uses  because 
(unlike  a)  we  have  a simple  explicit  expression  for  its  variance  (3.22) 
and  its  efficiency  is  competitive  for  a more  than  (say)  four  or  five. 
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where  i{i(a)  is  the  psl  function  (derivative  of  log  gamma).  The  maximum 

likelihood  estimators  (a ,8)  is  found  from  (4.2)  by  replacing  x with  x 
1 r n 

and  In  x with  In  x = — ) , In  x..  Thus  it  fits  into  our  generalized 

n ^1  X 

moment  scheme  utilizing  the  arithmetic  and  geometric  means.  To  solve  the 
system,  one  sets  3 = x/a  into  the  second  member  of  (4.2)  and  search  for 
the  root  of  In  a - ’>(a)  = In  x - In  x,  which  requires  a psi  function 
capability  [1,3].  This  is  not  difficult  with  a large  computer,  but  may 
be  a challenge  for  a small  one,  e.g.  a hand  held  calculator.  Viable 
alternatives  are  available  using  the  ordinary  method  of  moments  and  a 
• generalization  that  exploits  the  harmonic  mean. 

The  information  matrix  (2.9)  is 

/ a/8^  1/B  j 

A = (4.3) 

( 1/3  <i’(a)  ) 

whose  inverse  is 

\ 


aijj  ' (a  )-l 


- 6 


- 6 


(4.4) 


where  primes  denote  derivative. 

- 2 2 
The  ordinary  method  of  moments  equates  x and  s to  a6  and  a3  , 

respectively,  and  the  resulting  estimator  is 


6 = s^/x  , 


-2,  2 

a >•  X /s 


(4.5) 


Using  = 6 and  02  = a in  order  to  conform  with  (4.2),  one  can  apply 


(2.2)  to  obtain 


( “ 

( 2ot6 


e I 
) 


(4.6) 


Note:  Although  the  method  of  moments  shares  an  equation  with  maximum  likeli- 

hood  (x/B  - a/s  = 0),  there  is  no  advantage  to  using  this,  through  (2.18), 
in  this  case. 

The  matrix  C is  obtained  from  (C.5),  (C.6),  and  (C.7).  Then  (2.7) 
can  be  applied  to  get 


61  2a  + 3 
a 2(a+l) 


M = 2(a+l) 


-6 


- 6 


(4.7) 


and 


|mI  = 2(a+l)S^ 


(4.8) 


4 


I 


Turning  to  another  choice,  let  us  recognize  that  the  ordinary  method 
of  moments  uses  moments  of  order  one  and  two,  while  maximum  likelihood 
uses  moments  of  order  one  and  zero.  Heuristically  one  might  find  advantage 
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Ll 


in  moments  of  order  one  and  minus  one.  Consider  for  the  system  (2.2) 
(see  (C.2)  and  (C.8)), 


X - ae  =0 


-sir:!)  ■ ° 


where  y “ ^ I ^ 1/x^  and  the  consequent  estimators 


(4.9) 


* - 1 
6 = X - - 

y 


xy  - 1 


(4.10) 


* * 

which  satisfies  6 > 0 and  a > 0 since  the  arithmetic  mean  is  larger 
than  the  harmonic  mean.  Proceeding  as  before  we  obtain,  for  a > 2, 


i e^(a-] 


S(a-l)' 


(4.11) 


* _ 2(g-l)‘ 

a - 2 


2 I 2(a-l)^ 


2 2a  - 3 


(4.12) 


1m*1  = 28^ 

' ' a - 2 


for  a > 2 


(4.13) 


The  asymptotic  efficiencies  of  the  estimators  (4.5)  and  (4.10)  are 
computed  from  (2.22)  using 


|A|  = (am'M  - 1) 
3 


(4.14) 


These  efficiencies  do  not  depend  on  the  scale  parameter,  S,  and  can  be 


plotted  as  functions  of  a,  as  Is  done  in  Figure  4.1,  Also,  the  former 

* * --  - 

appears  in  [15 ]•  The  alternative  (a  ,8  ) catches  up  to  (a, 3)  at  a = 3 


and  remains  better  thereafter.  The  relative  efficiency  of  the  latter  with 
respect  to  the  former  is  (see  (4.13)  and  (4.8)) 


I * . 2 

|M  I , (g-1)^ 

iMl  (ct-2)(ct+l) 


(4.15) 


It  is  less  than  100%  for  all  a > 3 and  falls  below  90%  for  4 < a < 7. 


Figure  4.1 

Asymptotic  Efficiency 
( Gamma  ) 


Eff  (»,  yS  ) 


rs/  rsj 


Eff  (oi,  /3  ) 


6 7 8 9 10 


FIGURE  4.1.  Asymptotic  Efficiency  (Gamma) 


As  a numerical  example,  the  three  estimators  were  applied  to  135 
observations  on  the  service  time  of  the  check-in  queue  at  a major  hotel. 
The  gamma  distribution  fits  quite  acceptably  according  to  chi  square 
criteria.  The  results  are: 

I = .303  B = -305  3*  = -298 

a = 6.33  ot  = 6.31  a = 6.45 

Counterexample : This  opportunity  is  taken  to  show  that  the  representation 

(2.18),  which  is  the  main  theorem  of  [12]  , is  false  when  one  drops  the 
assumption  that  the  subsystems  u = 0 (2.11)  are  likelihood  equations. 

We  use  the  system 

y -?(^=  ° 

s"  - aS^  = 0 


which  has  an  equation  in  common  with  the  method  of  moments  (4.5),  u . this 
common  equation  is  not  part  of  the  maximum  likelihood  system.  This  system 
has  explicit  solutions 

2 I 


i 


+ 4s  / , 

y ) 


'*  1 
a = 1 + — 

y6 


One  readily  develops,  using  (4.11)  and  (4.6) 


S^(a-l) 


!(c*-l)‘ 


2aS 


and  from  (C.7),  (C.9),  (C.ll) 
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.>»  • * 


0 


3 


e^(a-l)^  (a-2) 


2a3  (a+3) 


Using  M ^ = A'C  ^A,  we  calculate 


6^(a-l)^(a-2)  8^(a+3) 


g-2  + 1 

8(a-l)  S(a+3) 


e(a-l)  B(a+3) 


■ iL-2„  + ^ 

(a-l)2  2g(a+3) 


and  this,  clearly,  has  no  submatrix  in  common  with  (upon  inverting  (4.7)) 


2a  + 3 
a 2(a+l) 


Negative  Binomial 


The  negative  binomial  density  is  parametrized 


, r(r+x)  X r 
f(x;r,p)  = q P 


£or  X = 0,1,2,...  (4.16) 

0 < r , 0 < p < 1 , p+q  = 1 


and  the  partial  derivatives  of  its  logarithm  are 


Sp  = r/p  - x/q 


= ii(r+x)  - ;(r)  + In  p 


Using  the  basic  recursive  formula  for  the  psi  function  [1], 


(4.17) 


.'3* 


■ • ■ 


X ^ 

^(r+x)  - ^(r)  = I - -_j 

j=l  ^ 


one  may  express  the  system  of  maximum  likelihood  equations  as 


X - rq/p  = 0 


(4.18) 


1 

\ r---r+- 

1=1,..., n j=l  ■’ 


+ In  p = 0 


Solving  (4.18)  is  quite  difficult  to  manage  without  a large  computer. 
Appendix  F contains  an  APL  program  (PSIB)  which  accomplishes  it. 

The  information  matrix  can  be  developed  using  (2.9).  The  result  is, 
using  6^  = p and  62  = r. 


I -1/p 


(4.19) 


where 


A 22  = '/(r)  - E(^’(r+X)) 


(4.20) 


The  properties  of  A22  are  developed  in  Appendix  D (see  (D.22)  and  (D.26)), 
along  with  the  series  representation  of  the  determinant  (D.25), 


, ” n 

= Y _9 

2 4 u + 1 
p n=l 


(n+r) : 


(4.21) 


which  converges  rapidly  for  most  values  of  p.  For  purposes  of  computation, 
one  may  as  well  use  (4.21)  in  (4.19)  and  solve  for  ^2'}'  thus 


*22 


(4.23) 


To  estimate  by  the  ordinary  method  of  moments,  we  use  (D.5)  and  form 


Let  us  turn  to  the  question  of  using  a different  moment  in  the  second 


equation.  Since  the  negative  binomial  has  positive  probability  mass  at 

zero,  we  use  averages  of  the  (1+x^)  ^ as  was  done  when  the  Poisson  case 

was  considered.  Perhaps  that  level  of  success  can  be  matched. 

Letting  y = — Z,  l/(l+x,),  use  for  the  system  (2.1) 
n J.  1 


px  - rq  = 0 

(r-l)qy  - p + p^  = 0 


(4.30) 


because  of  (D.7).  The  system  (4.30)  cannot  be  solved  explicitly,  but  it 
can  be  managed  with  a hand  held  calculator.  Use  the  first  member  to  obtain  r 
as  a function  of  p and  its  derivative 


dr  _ X 
dp  2 

q 


(4.31) 


and  substitute  into  the  second  member  of  (4.30)  to  obtain  a function  f 
of  the  form 

f(p)  = p*^  + (r-l)qy  - p 


= p - y + p (y  - 1 + xy) 


(4.32) 


and  having  derivatives 


a 

dp 


(y  - 1 + xy) 


p^x(q  + In  p) 

q" 


* 


I 


4 


The  solution  of  f(p)  = 0 can  be  obtained  by  Newton's  method,  always 
remembering  to  update  r as  well  as  p.  The  initialization  p = .5  and 
r = X appears  to  be  satisfactory,  but  normally  convergence  would  be  faster 
if  the  moment  estimators  (4.25)  are  used. 


22 


trtri'i>iTA-|i 


% 

=« 


1e  ic 

The  resulting  estimator  will  be  denoted  by  p , r . From  (4.32)  it 
is  seen  that 

f(0)  = -y  < 0 , f(l)  = xy  > 0 

•k  k 

so  that  0 < p <1,  and  r >0  follows  from  this  by  (4.31). 

Let  us  turn  to  the  development  of  the  asymptotic  covariance  structure 

* k 

of  p , r . Taking  partial  derivatives  and  the  expectation  of  (4.32)  to 
produce  the  coefficient  matrix  A of  (2.2)  yields 


-q 


In  p 


1 

1 


(4.33) 


with  the  help  of  (D.7).  Attention  is  drawn  to  the  fact  that 


limit  ^ ~ P = p - p In  p 
r -*■  1 ^ 


(4.34) 


The  covariance  matrix  C of  (2.6)  is  derived  from  (4.30)  with  the  help  of 
(D.5),  (D.7),  and  (A. 5).  The  result  is 


rq 


r /-I 

rqp  - p(l-p  ) 


\ . 

'rqp  - 


rqp*^  - p(l-p^) 


>1 


(4.34) 


Let  us  draw  attention  to  the  fact  the  first  equation  of  (4.17)  is  a multiple 
of  the  first  equation  of  (4.32)  rather  than  being  identical.  This  is  the 
reason  that  (4.33)  and  (4.35)  are  modifications  of  (2.17)  and  (2.16)  rather 
than  exact  analogies.  Thus  the  bookkeeping  that  follows  must  be  done 
carefully. 
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The  direct  development  of  the  matrix  M from  (2.7)  with  (4.34)  and 
(4.35)  as  input  is  messy.  Instead  let  us  recognize  that  (4.30)  shares  an 
equation  with  the  likelihood  system  (4.18)  and  use  (2.18).  Thus,  with  A 
given  by  (4.19),  we  have 


(4.36) 


r/qp^ 

-1/p 

-1/p 

22 

M 

where  is  obtained  from  (2.23)  and  |c|  from  (2.16). 


I r 2 2 

I 2 ®22  p ®21®22 

[qp 


(4.37) 


1^1  “ 2 ^22  ~ ®21 

qp 

and  (§21’  ®22^  second  row  of  A in  (4.33)  and  is  taken  from 

(4.35).  Using  (2.24)  we  obtain 


1*1-1  _/  r ^ 1 
|M  I - 1 2 ®22  p 

V qp 


(4.36) 


|c|  = ^ (r-1)^  - 4l 


The  asymptotic  efficiencies  of  (p,r)  and  (p  ,r  ) appear  in  Tables 

4.1  and  4.2,  resp.,  for  p = .1(.1).9  and  r = .5(.5)5,  6(1)19  where  the 

parentheses  indicate  the  Indices  of  advancement.  The  efficiency  of  (p,r) 

is  monotone  increasing  in  r for  each  p.  It  is  lower  for  the  smaller 

* * 

values  of  p.  The  efficiency  of  (p  ,r  ) is  high  for  the  smaller  values 
of  r and  decreases  (generally).  It  is  not  monotone  for  p = .1,  .2. 


* * ^ - 

The  relative  efficiency  of  p , r with  respect  to  p,  r,  i.e., 

Rel  eff  = 

A “ff 

appears  in  Table  4.3.  Generally  p , r is  preferable  for  r less  than 
or  equal  to  (say)  2.5.  In  general  p,  r is  preferable  elsewhere,  but  it 
does  not  matter  much  for  small  values  of  p. 

The  three  estimation  schemes  were  applied  to  the  Cricket  score  data 
of  Reep,  Pollard  and  Benjamin  [13],  which  provided  the  following  comparisons. 


Cowdry 

Barrington 

Graveney 

X 

1.692 

2.095 

1.570 

2 

s 

4.343 

4.939 

4.474 

(l+x)"^ 

.603 

.538 

.626 

n 

156. 

116. 

107. 

P 

.329 

.345 

.317 

r 

.831 

1.014 

.729 

P 

.390 

.424 

.351 

r 

1.080 

1.543 

.849 

* 

P 

.371 

.326 

.389 

* 

r 

1.000 

1.012 

1.000 
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Beta 


The  beta  density  has  the  form 


f(x;a,6)  = r(a“)'^r(e)  ^(1-x)® 


(4.37) 


for  0 < X < 1,  0<  a,  0 < 6 


and  the  partial  derivatives  of  its  logarithm  are 


= iJ;(a+B)  - ip(a)  + In  x 
S = i|;(a+B)  - + In(l-x) 

P 


(4.38) 


The  information  matrix  is,  for  0^^  = a,  = B 


4)'(a)  - i|;'(oH-B) 


- t|)(cri-B) 


-ijj'  ((^B) 


4*' (6)  - i>'(ofS) 


(4.39) 


The  system  of  maximum  likelihood  equations 


In  X = i)'(a)  - 4'(ori-6) 


In(l-x)  = i(j(6)  - 4/(a+6) 


(4.40) 


uses  the  geometric  means  of  x and  1-x,  and  is  difficult  to  solve.  UTiat 
other  pairs  of  statistics  might  be  substituted? 

Clearly  x and  (l~x)  cannot  be  paired  since  they  are  functionally 
related.  The  latter  is  merely  1-x.  Let  us  first  develop  the  ordinary 
method  of  moments.  Using  (E.4)  and  (E.5),  choose  the  systems 


(a+e)x  - a = 0 


(a+6)^  (a+e+l)s^  - aB  = 0 


(4.41) 


and  solve  explicitly  for 


'..i«0KV 


V -4-f 


^ - S - x(l-x) 


s - x(l-x) 


(4.42) 


It  may  occur  that  ot  < 0,  6 < 0. 

The  coefficient  matrix  A of  (2.2)  is 


(4.43) 


a + 6 + 1 


(4.44) 


and  C of  (2.6)  takes  the  form 


a8 

a+6+1 


(a+B)  "^(a+B+l)  Cov(x 


(a+S) ^(a+S+1)  Cov(x,s^)  (a+6) ^ (a+6+1) ^ var 


x,s^)  I 

(s^)  j 


(4.45) 


and  the  use  of  (E.2)  thru  (A.l)  and  (A. 2)  does  not  appear  to  simplify  in 

II  I I I 1 2 

any  useful  way.  Appendix  F contains  programs  to  compute  |M|  = |C|/|A|‘'. 
Let  us  turn  to  the  pair  of  statistics 


1 f 1 

^ ° n ^ r 
n 1 Xi 


-I  — 

n { 1-x. 


(4.46) 


which  are  not  functionally  related.  Using  (E.8)  we  may  choose  the  system 


(a-l)y  - (a+3-1)  = 0 
(3-1) z - (a+S-1)  = 0 


(4.47) 


for  a 1,  3 > 1.  The  solution  is 


(4.48) 


« 


I 


= y(z-i) 


z(y-i) 


yz  - y - z 


yz  - y - z 


Clearly  y > 1 and  z > 1 but  the  denominators  are  not  necessarily  positive 
since  zy-y-z+l=  (y-l)(z-l)  can  be  less  than  one. 

For  the  system  (4.47)  the  coefficient  matrix  becomes 


( A- 

-1 

* 1 

A =; 

1 -1 

a 

V 

6-1 

and  (E.IO) 

one  can 

(4.49) 


(a+2-1) 


(4.50) 


for  a -■  2,  3 > 2.  Use  of  (4.49)  and  (4.50)  in  (2.7)  does  not  produce  a con- 

* 


venient  expression  for  M . However  its  determinant  is  easily  managed. 

For  fun,  let  us  also  try  a system  based  on  moments  of  order  one  and 
minus  one.  CThis  would  seem  to  straddle  the  two  geometric  means  appearing 
in  maximum  likelihood.)  We  are  lead  to 


(a+B)x  - a = 0 
(a-l)y  - (a+6-1)  = 0 


(4.51) 


and  the  resulting  estimate 


-*  (v-l)x 
a = — “ 


yx-1 


(y-D(l-x) 

O — 

v^-1 


(4.52) 


1 


and  this  satisfies  a > 0,  6 >0.  (But  do  not  forget  that  use  of  (E.6) 

in  (4.51)  requires  a > 1.)  Proceeding  in  the  usual  way,  we  calculate 


6 

a + B 


aB 

a + B + 1 


a + B 


B(g+6+l) 
a - 2 


(4.53) 


4.54) 


Again  M does  not  have  a convenient  form.  The  determinates  of  (4.52)  and 
(4.53)  are  expressed 


-3 

(a+B) (a-1) 

2B^ 
a - 2 


(4.55) 


for  a > 2. 

The  efficiencies  of  (4.42),  (4.43),  and  (4.52)  are  compared  in  the 
tables  that  follow.  Table  4.4  contains  the  efficiency  of  the  ordinary  moment 
estimator.  Efficiency  is  high  if  a is  not  too  far  from  3 and  both  are 
at  least  two.  Elsewhere  they  are  low,  but  still  the  choice  because  all  the 
numbers  in  Table  4.4  are  better  than  their  competitors  in  Tables  4.5  and  4.6. 
The  pair  (a  ,S  ) is  generally  better  than  (g  ,B  ) but  not  uniformly  so. 

It  matters  little  since  (a, 3)  is  the  "hands  down"  winner.  This  result 
parallels  what  was  learned  in  the  symmetric  beta  case.  The  variable  X ^ 
is  unstable  in  this  population. 


TABLE  4.4.  Eff(a.3)  BETA  POPUI.ATION 


.5 

1.0 

1.5 

2.0 

2.5 

3.0 

3.5 

.5 

0-.  4 93 

0.554 

0.537 

0.507 

0.477 

0.451 

0.42  9 

1.0 

0.554 

0.713 

0-747 

C.739 

0,719 

0.696 

0.674 

1.5 

0. 537 

0,  747 

0.820 

0. 839 

0.  835 

0.822 

0.8  06 

2.0 

0.507 

0.739 

0.839 

0.878 

0.3  89 

0.367 

0.378 

2.5 

0 . 477 

0.  71  9 

0.835 

0.889 

C.912 

0.919 

0,918 

3.0 

0 . 451 

0. 696 

0.822 

0.887 

0.919 

0.934 

0.938 

3.5 

0.429 

0 .674 

0.806 

0.878 

0 . 91  8 

0.938 

0.948 

4.0 

0.  41 1 

0.  653 

0.790 

0.  867 

C.912 

0.938 

0.952 

4.5 

0.  395 

0.  635 

0.774 

0.855 

0.904 

0.933 

0.951 

5,0 

0.  382 

0.  618 

0.758 

0.843 

0.895 

0.923 

0.948 

5.5 

0.  370 

0.  604 

0.744 

CO 

o 

0.886 

0 , 921 

0.944 

6.0 

0.  360 

0.591 

0.731 

0. 81  9 

0.876 

0,914 

0.938 

6.5 

0.351 

0.579 

0.719 

0.809 

0.867 

0.906 

0.933 

7.0 

0. 344 

0.568 

0.709 

0.799 

0.858 

0.8  99 

0.927 

4.0 

4.5 

5.0 

5.5 

6.0 

6.5 

7.0 

.5 

0.411 

0.395 

0.382 

0.370 

0.360 

0.  351 

0.344 

1.0 

0.653 

0.635 

0.618 

0.604 

0.591 

0.579 

0.568 

1.5 

0.790 

0.774 

0.758 

0.744 

0.731 

0.719 

0.709 

2.0 

0.867 

0.855 

0.843 

0.831 

0.819 

0.809 

0.  799 

2.5 

0.  912 

0.904 

0.895 

0,666 

0.876 

0.  367 

0.858 

3.0 

0 . 93  8 

0.  933 

0.928 

0.  921 

0.914 

0.  906 

0.899 

3.5 

0 .952 

0.  951 

0.943 

0.944 

0.  93  8 

0.  933 

0.927 

4.0 

0.959 

0.961 

0.961 

0.  95  6 

0.955 

0.951 

0.946 

4.5 

0.961 

0.966 

0.  968 

0.968 

0.  966 

0.  963 

0.960 

5.0 

0 . 961 

0 . 968 

0.972 

0.  973 

0.  973 

0.972 

0.969 

5.5 

0. 958 

0, 968 

0.973 

0.976 

0.  977 

0.  977 

0.  976 

6,0 

0 .955 

0.966 

0.  973 

0.977 

0.980 

0.  980 

0.  980 

6,5 

0 . 951 

0.963 

0,972 

0.  977 

0.980 

0.  982 

0.  983 

7.0 

0.946 

0.960 

0.969 

0.976 

0.980 

0.  983 

0.  934 

4 


I 


33 


V 


O' 

ro 

<£> 

to 

iO 

C7> 

CO 

o 

<o 

CM 

r^ 

o 

CM 

fO 

tn 

in 

r^ 

tn 

r- 

00 

OD 

o> 

CT' 

cr> 

(T> 

cr> 

<D 

o 

o 

o 

o 

o 

O 

o 

O 

o 

O 

O 

a* 

t/> 

tn 

m 

ro 

to 

3- 

o 

:3* 

fO 

CM 

o 

CM 

r> 

=r 

LO 

lO 

m 

00 

00 

Oi 

0> 

a> 

cr> 

CT> 

0^ 

vO 

• 

* 

• 

• 

• 

• 

• 

• 

• 

• 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

O 

CO 

csj 

o 

C7> 

3* 

CO 

CM 

o 

rH 

CO 

CO 

J- 

cf) 

CO 

CO 

o 

0> 

o> 

o> 

<n 

CJ> 

sO 

O 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

CM 

fH 

o> 

o 

CO 

cr 

rH 

cO 

a> 

ro 

CM 

CO 

<n 

CM 

CO 

ro 

CO 

• 

m 

r' 

00 

CD 

00 

(D 

(J> 

(J) 

at 

03 

lA 

* 

« 

• 

• 

• 

• 

o 

o 

o 

O 

o 

o 

o 

O 

o 

o 

o 

M 

H 

<J 

2 

o> 

<Ji 

in 

ro 

ro 

cn 

ro 

03 

CO 

CM 

lO 

CO 

O 

rH 

rH 

(M 

CM 

O 

tn 

00 

00 

00 

o> 

CJ» 

CO 

C73 

03 

PU 

lA 

« 

• 

• 

• 

♦ 

• 

• 

• 

• 

• 

o 

o 

o 

o 

o 

o 

O 

o 

O 

o 

E- 

;sj 

CO 

m 

in 

CO 

CO 

t^3 

CM 

in 

ro 

CM 

o 

cr 

CO 

O 

o 

o 

o 

tn 

CO 

CO 

03 

00 

CO 

03 

03 

03 

ca 

o 

o 

o 

o 

o 

o 

O 

o 

o 

o 

‘ « 

cD 

o 

to 

CD 

rH 

<T> 

CO 

tn 

u 

c 

CM 

O 

cc 

(N 

4- 

O 

p*- 

p' 

in 

CO 

CO 

CO 

CO 

CO 

CD 

CO 

vr 

• 

• 

♦ 

• 

« 

• 

• 

• 

• 

* 

lO 

o 

o 

o 

O 

o 

o 

o 

o 

o 

o 

<r 

UJ 

to 

o 

<T» 

to 

tn 

*n 

rH 

4* 

n 

uO 

o 

CD 

o 

CM 

CM 

CM 

CM 

< 

m 

tc 

C*^ 

00 

CO 

CO 

CO 

CO 

CO 

H 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

o 

o 

o 

o 

o 

o 

o 

o 

o 

CM 

CO 

to 

to 

rH 

<73 

CM 

4' 

CO 

o 

to 

r^ 

o 

CM 

CM 

ro 

CO 

CO 

CO 

• 

4“ 

tO 

to 

p^ 

P^ 

p** 

p** 

p^ 

<r» 

• 

* 

• 

• 

* 

• 

• 

o 

CO 

o 

o 

o 

o 

o 

o 

o 

o 

r 


TABLE  4.6.  Eff(a  ,S  ) BETA  POPULATION 


.5 

1.0 

1.5 

2.0 

2.5 

3.0 

3.5 

2.5 

0.119 

0.  202 

0.258 

0.297 

0.326 

0.347 

0.364 

3.0 

0 .1  63 

0 . 276 

0.  353 

0.407 

0.447 

0.477 

0,5  00 

3 . 5 

0.184 

0.313 

0.400 

0.461 

0.506 

0.540 

0.566 

4.0 

0.196 

0.  333 

0.426 

0.492 

0.540 

0.  576 

0 . 605 

4.5 

0.204 

0.346 

0.  443 

0.511 

0.561 

0,599 

0,629 

5.0 

0,209 

0.355 

0.454 

0,524 

0.576 

0.615 

0.645 

0.  21  2 

0.  361 

0. 462 

0.534 

0.586 

0,626 

0.657 

5 . 5 
6.0 

0 . 21  5 

0.  366 

0.  468 

0.541 

0.594 

0.634 

0.666 

0.  217 

0.369 

0.473 

0.546 

0.  600 

0.641 

0.673 

O . J 

7.0 

0.  21  8 

0.372 

0.476 

0.550 

0.604 

0,  645 

0.678 

4.0 

4.5 

5.0 

5.5 

6.0 

6.5 

7.C 

2 , 5 

0.378 

0.  389 

0.398 

0.406 

0.412 

0.  41  8 

0.423 

3.0 

0.519 

0.534 

0.546 

0.557 

0.566 

0.574 

0,581 

3.5 

0,588 

0. 605 

0. 61  3 

0.632 

0.642 

0.  651 

0.659 

4.0 

0,627 

0.646 

0.662 

0.  675 

0.686 

0.  695 

0.704 

4.5 

0.653 

0.672 

0.688 

0.702 

0.  714 

0.  724 

C.  735 

n 

0.670 

0.  690 

0.707 

0.721 

0.733 

0.743 

0.752 

5.5 

0.682 

0.  703 

0.720 

0.734 

0.746 

0.757 

0.766 

6.0 

0.691 

0.712 

0.730 

0.744 

0.757 

0.768 

0.777 

6. 5 

0.698 

0.  719 

0.  737 

0.752 

0.765 

0.776 

0.785 

7.0 

0.704 

0.725 

0.743 

0.75  8 

0.771 

0.782 

0.792 

. • . 


(A. 2) 


T 


« 

=» 


2 1 2 4 2 4 

Var(s  ) = — {ffl,  - 4m_mT  + 3m_  - 4o  + — r o } 
n 4 j i z n'^i 


1 , 4,2  4, 

— lu.  - O + :r  a } 

n 4 n-1 


2 


(u^-2u2> 

2 

n 


+ 


(A. 3) 


Cov(X.Y)  = - {1  - E(X)  E(l/X)}  < 0 
n 

Cov(X.Y’)  = ^ {1  - [1  + E(X)]  E(^)}  < 0 

Cov(s^,Y)  = - - {m,E(l/X)  + m - 2m?E(l/X)} 
n z 11 

= - , - u^m  , + u) 

n Z -1  -1 

Cov(s^.Y')  = - J {[^2  - ^ 


(A.4) 
(A.5) 
(A. 6) 

(A. 7) 


The  proofs  of  (A.l)  to  (A. 7)  follow.  It  is  convenient  to  record 
the  relationships 

2 

“2  “ *^2 

+ 3u2iJ  + (A. 8) 

2 4 

+ 4iijU  + 6^2^  + U 

To  enhance  the  readability,  the  symbols  E,  V,  C will  be  used  to  denote 
expectation,  variance,  and  covariance  (resp.).  Parentheses  and  subscripts 
will  be  used  sparingly. 
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Proof  of  (A. 4).  Consider 


m EEX^Ed/X  ) 

- ^ {n  + n(n-l)  EXE(1/X)} 
n 


from  which  we  subtract  EXE(1/X)  to  produce  (A. 4).  The  fact  that  (A. 4)  is 
negative  follows  from  the  presumption  that  X > 0 and  the  consequence  that 


the  harmonic  mean  is  less  than  the  arithmetic  mean.  Thus  (ave  X. ) x (ave  > 1 

1 

unless  (all  = constant)  and  apply  the  law  of  large  numbers. 


Proof  of  (A. 5).  Follows  from  the  fact  that  C(1  + X)Y'  = CXY'  and  the 
application  of  (A. 4). 


Proof  of  (A. 6).  Consider 


n(n-l)  Es^Y  = E E Z (X.-X)^ 

Xf  3 

1 2-2 

= E Z (Exf-nX^) 

\ 3 

X^  . X.X. 

= EEE  -1  - i EEE 

Xf  n X^ 


= nEX  + n(n-l)  EX^E  - EX  - 2(n-l)  EX 

- (n-1)  EX^E  I - (n-1) (n-2)  E^XE  ^ 


and  divide  through  by  (n-1)  to  get 
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E(z:X^)^  = nEX^  + n(n-l)E^X^ 


(A. 9) 


EX^(ZXj)^  - nEX^  + n(n-l)EV  + n(n-l)  (n-2)EX^E^X  + 2n(n-l)EX^EX  (A. 10) 


ECrx^)^  = nEX^  + 4n(n-l)EX^EX  + 3n(n-l)E^X^  + 6n(n-l) (n-2)EX^E^X 

+ n(n-l)(n-2)(n-3)E'^X  (A.  11) 


The  proper  combination  of  (A. 9),  (A. 10),  and  (A. 11)  produces 


Es^  =“  - {EX^  - AEX^EX 


+ [(n^-2+3)E^X^-2(n-2)(n-3)EX^E^X+  (n-2)  (n-3)E^x] } 
n~i 


The  subtraction  of  o in  the  form 


4 2 2 2 2 4 

o - E X'^  - 2EX'‘e  X + E X 


yields 


2 4 'K  2 2 424 

nVs  = EX  - 4EX^EX  + 3E  X - 4o  + o 

n-1 


and  this  is  (A. 2).  The  second  version  follows  from  using  (A. 8)  and  modifying. 

Proof  of  (A. 3).  The  variance  of  this  form  of  the  sample  variance  can  be 
developed  from  (A. 2)  using  (A. 8).  It  also  appears  in  [7,  p.  183]. 


APPENDIX  B 


Moments  of  the  Poisson 


The  Poisson  random  variable  X had  density 


f(x;X)  = e ^ X^/xI  , X = 0,1,2,. 


(B.l) 


and  it  is  well  known  that  u = A,  o = X.  The  probabi] '.ty  generating  function 
is 

G(u)  « E(u^)  = (B.2) 


I 


and  the  moment  generating  function  can  be  obtained  from  (B.2)  by  the  replace- 
ment e'^  for  u.  Moments  can  be  obtained  by  repeated  differentiation. 

We  record 

E(X)  = X 


E(X^)  = X^  + X 


E(X^)  = X + 3X^  + X^ 


4 2 3 3 

E(X  ) = \ + 7X  + 6X-^  + X-^ 


Use  of  (3.3)  into  (A. 2)  produces 


= X + 3X' 


Var(s^)  = - {2X^  + X}  + o(-) 
n n 


Cov(X,s^)  = — X 
n 
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V-** 


(B.3) 


(B.4) 


(B.5) 


(B.6) 


The  moments  of  (1+X)  ^ can  be  obtained  by  Integrating  the 
generating  function  (B.2) 


E(t^)  “ E / u^du  = / G(u)  du  * e e^“  ‘Ju  “ f (B.7) 


E(-r^)^  E / / du  dv  = ff  G(uv)  du  dv  = e ^ du  dv 

^ * 0 0 

-X  1 Xu  , -X  “ 1 j-1 

= V / ^ ‘i-  V I / ^ 


^ b " 


^10  ^ 


-X 


j-1 


y ^ 

L 4 . i I 

1 J 


(B.8) 


This  opportunity  is  taken  to  record  an  alternative  way  of  obtaining 
moments,  (B.2),  which  in  this  case  is  somewhat  easier  than  the  differentiation 
of  the  moment  generating  function.  One  begins  with  the  generating  function 
G(u)  = E(u  ) and  replaces  the  argument  u by  a product  of  dummy  variable 
uv  • • • z containing  as  many  factors  as  the  order  of  the  moment  to  be  calcu- 
lated. Then  one  takes  a partial  derivative  with  respect  to  each  variable 

u,  V,  etc.  and  the  desired  moment  is  obtained  when  each  variable  is  set  to 

2 

unity.  For  example,  E(X  ) can  be  obtained  from 


uv  3u  3v 


G(uv) 


3u  3v 


E(uv)^  = E(X^u^“^v^“^)  (B.9) 


The  four  moments  (B.3)  can  be  obtained  in  this  way  replacing  u 
with  uvwz  in  (B.2).  The  resulting  derivatives  are 
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•**  A' 


G » AG 
u 


Ik. 


= XG[1+  Auvwz] 

2 2 
G = XG{z[l  + Xuvwz]  + Xuvwz  [1  + Xuvwz]  + Xuvwz  } 
uvw 

G = Xuvw  G + XGffl  + Xuvwz]  + zXuvw  + 2zXuvw(l+Xuvwz) 
uvwz  uvw 

^,2  2222_^_,  I 

+ X u V w 2 + 2zXuvwl 


J 


APPENDIX  C 


Moments  of  the  Gamma 


The  gamma  random  variable  X has  density 


f(x;a,6)  = ^ ^ x“ 

r(a)e“ 


(C.l) 


and  It  Is  well-known  that 


U = a6 


2 «2 
o = ag 


(C.2) 


Direct  Integration  produces  the  formula,  for  r > 0 


p/yf-i  = r (g+r) 

E(X  ) - 8 


(C.3) 


and  for  r < g 


^ r(g) 


(C.A) 


Use  of  (C.3)  In  (A.l)  and  (A. 2)  produces 


Var(X)  = - gB^ 
n 


(C.5) 


- 2 2 3 

Cov(X,s  ) = — g8 
n 


(C.6) 


Var(s^)  = 2g6^ 

n-1  n 


(C.7) 


« 


Using  (C.A)  one  readily  calculates 


I 


AA 


, -rf  jM,  ■» 


• • ■ 


rr  r 


E(^)  = ^ 


'X'  8(a-l) 


1 < a 


(C.8) 


6 {a-l)"(a-2) 


2 < a 


(C.9) 


and  letting  ^ ^ I”  (1/X^),  we  find  that,  with  the  aid  of  (A. 4), 


Cov(X,Y)  = 


-1 


n(a-l) 


1 < a 


(C.IO) 


and  with  the  aid  of  (A. 6) 


Cov(s  ,Y)  = 0 


(C.ll) 


Because  of  the  formula 


r'(y)  = / ln(y)  y“  ^ e ^ dy 
0 


(C.12) 


one  can  develop,  for  r > - a, 


E(X^  In  X)  = (In  B + iKa+r)) 


(C.13) 


since  r'(oc)  = r(oc)  'Jj(a), 


4 


I 
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APPENDIX  D 


Moments  of  the  Negative  Binomial 
The  density  has  the  form 

f(x;,r,pj^  x'r(r)  ^ ^ for  x — 0,1,.. • (D.l) 

0<r,0<p<l,  p+q  = 1 

Its  probability  generating  function 

Y r 

G(u)  = E(u^)  = (D.2) 

(1-qu)*^ 

will  be  exploited  broadly.  Successive  derivatives  of  (D.2)  evaluated  at 
u = 1 produce  the  factorial  moments 

= EX(X-l)  •••  (X-s+1)  = (q/p)®  r(r+l)  •••  (r+s-1)  (D.3) 

to  which  one  may  apply  some  orderly  substitution  and  obtain  the  first  four 
moments,  using  A = q/p 

E.\  * Ar 

rx~  - A^r(r+1)  + Ar 

(D.4) 

EX^  - A^r(r+l)(r+2)  + 3A^r(r+l)  + Ar 

E.i'^  • -4‘^r(r+l)(r+2)(r+3)  + 5A^r  (r+1)  (r+2)  + 4A^r(r+l)  + Ar 
The  mean  anci  variance  of  che  population  are 


^6 


V = rq/p 


(D.5) 


2 , 2 
o = rq/p 


Use  of  (D.4)  in  (A.l)  and  (A. 2)  provides  the  covariance  matrix  for 
the  ordinary  method  of  moments 

n Var(X)  = A^r  + Ar  = rq/p^  = o“ 


n Cov(X,s^)  = 2A^r  + 3A^r  + Ar  = 

P 

n Var(s^)  = 2A^r(r+3)  + 4A^r(r+3)  + A^r(2r+7)  + Ar 

= o^[2(r+3)qVl/p^ 


(D.6) 


The  harmonic  mean  alternative  requires 


U+X''  (r-l)q 


JL)2  . e‘  i r k_  . i] 


for  r ^ L.  The  case  r = 1 is  treated  in  [12,  Appendix  A].  Expressions 
(D.7)  and  (D.8)  are  justified  next  along  with  computational  formulae  for 
(D.8)  when  r is  either  a whole  number  or  a whole  number  plus  0.5. 


Proof  of  (D. 7) . The  generating  function  (D.2)  is  integrated. 

^^I+X^  ° ^ G(u)du  = p’^/  — 2 = ^ 

0 0 (l-qu)”'  (r-l)(-q)(l-qu)'' 


(r-l)q 


as  required. 


<4 


Proof  of  (D.8).  Replace  u by  uv  and  integrate  twice. 

- I I <=<”)  dv  . / d,  ] J 


. pI-  i f 1 l] 

(r-l)q  J,  uL  r-1  J 


as  required. 


Computational  formulae  require  managing  integrals  of  the  type 


j = r i 

0 Ud-qu)’^  “-1 


Clearly  = 0.  It  is  useful  to  apply  the  partial  fraction  representation 
repeatedly. 


uCl-qu)*^  (1-qu)^  u(l-qu)’ 


(D.IO) 


Since 


I .q-du  , JL  Lj^  . i] 

J ,r  r-1  r-1  -"J 

0 (1-qu)  i-n  -1 


we  have,  from  (D.IO), 


I = J + r — - ii 

r r-1  r-1  I r-1  J 


Let  <r>  be  the  greatest  integer  in  r and  apply  the  above  formula  to  get 


the  representation 


r r-j  r-1  r-< 

j*l  Lp  J J 


(D.ll) 


which  is  valid  for  r ><r>. 
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Let  r = < r>. 


One  can  show  directly  from  (D.9)  that 


, f q du 
J,  = J — r = - In  p 

1 (1-qu) 


(D.12) 


and  then 


\ 7^  [-jq  - i]  - i"  P - ^ 's'  J IP^'^-P'!  - IP  P ("-13) 

j=l  -J  Lp  J J p j=i  J 

o 

Accordingly  it  is  recognized  from  (D.8)  and  (J.9)  that  E(l/1+X)  = [p’"/ 'r-l)q] 

which,  together  with  (D.7)  and  (D.13)  enables 


''“'ife)  ■ J - p'  i"  p - I <“•“) 


for  r an  integer  ^ 2 (empty  sura  is  zero),  and  for  r = 1, 


Var(^)  = J { I 4 ■ J 

i+A  Cl  . - .2  q 

^ j-i  J ^ 


(D.15) 


This  latter  formula  is  developed  in  [12],  see  (A. 3)  and  (A. 5). 


Let  r = <r>  + .5.  The  exploitation  of  (D.ll)  requires  dealing  with 


'^r-<r>  '^1/2 


= /'[■  — ^ M du  . 

0 I-  1-qu  -* 


2 

Making  the  change  qu  = sin  9 provides  for  manageable  Integrals  of 
trigonometric  functions.  See  [14,  p.  316].  The  result  is 


T = T 1 n 
1/2  ^ 


1 + p p 


(D. 16^ 
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and 


E(^) 


<r> 


k"  p’'  hi 2^ 


1+X"  (r-l)q 


(D.17) 


which  is  valid  for  <r>  ^ 1.  For  r = 1/2  one  needs,  from  (D.9), 


"-1/2  ■ I 


/l-qu  _ 1 
u u 


du 


2 ^ 


qu 


u=i  1 r 1 

if  ^ 

-A] 

1 p 

u=0  0 *-u*'  l-qu 

^ J 

du 


2>/p  - 2 + J 


1/2 


(D.18) 


using  formula  [14,  p.316].  From  (D.8) 


^'1+X^  “ ~ q *^P  - 2 + J^,^}  = P-  In  ^7^1  (D-19) 

for  r = 1/2. 

The  information  matrix  (4.19)  contains  a difficult  element  •'22’ 
i-.lO).  It  may  be  managed  using  an  integral  representation  of  the  tri- 
gamna  function  [1,  p.  259], 

1 


. t / ^ r in  u r-i  , 

■i>  {r)  = - J Y^—  u du 

0 


(D.20) 


It  follows  that 


* 

4 


A22  = 'l-'(r)  - E{il;’(r+X)} 
1 


= - J " E(u^~^'^^)  ] du 

0 

= . / ^^r-1  r £l_  ^ 1 

0 L (l-qu)*^  -J 


(3.21) 


I 
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•4 


This  relationship  can  be  expressed  better  if  we  make  the  change  of  variable 


w = pu/(l-qu) 


and  manipulate  to  obtain 


ln(p+qw)  , r-1 


^22  J 1 - w 


(D.22) 


There  is  advantage  in  using  (D.22)  in  the  development  of  the  determinant  of 


A,  (4.19), 


1A|  ^ 

qp  0 P 


11.  i£(2±a«i  d(„") 

2 ^ q J 1 - w 


^ 

qp  0 


1 - w 


using  partial  integration  and 


ln(p+qw) 
1 - w 


= - ^ -3-  (1-w)’ 


In  this  form  it  is  easily  checked  that 


ii.it  !a|  . d. 

r -<■  0 qp  0 


-ln(p)~3, 

2 

qp 


limit  IaI  = 0 


For  computational  purposes  one  can  exploit  the  expansion 
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» n 


M2±awl)  . j al  (n-l)(l-w)"-2 


1 - w 


n»2 


which,  when  used  in  (D.23),  produces 


1 r n n-1  f r,,  » 

— o I q — J w 

^ H /\ 


n-2 


2 

qp  2 


dw 


1 r n n 
')  L a _ . 1 


P 1 


n n rl (n-1) I 
n+1  (n+r) J 


(D.25) 


1 “II 

1 r q rl  n! 

2 n+1  (n+r)I 

p n=l 


Similar  efforts  applied  to  (D.22)  can  produce 


22  ii„2(r4.-l): 


(D.26) 


APPENDIX  E 


Moments  of  the  Beta  Distribution 


The  Beta  random  variable  X has  density 


for  0<x<l,0<a,0<6 


r (a)  = / ^e  ^ dx 


This  is  called  the  B(a,6)  distribution  and  it  is  useful  to  note  that  1-X 
has  a B(6,a)  distribution. 

Moments  are  obtained  directly  by  manipulating  gamma  and  beta  functions. 


For  r > 0 


p/yr.  _ r(g-B)  r(g+r)  ^ (g4B-l).’  (g+r-1) 

^ “ r(g)  r(gH^+r)  (g-DI  (g-f«+  r-1)  I 


r r..  ,s,  _ r (g->6)  r (g+r)  F (3+s) 

t,lA  j - r(3)  r(g46+r+s) 


and  the  mean  and  variance  follow 


(E.4) 


g + 6 


(g+6)  (g-ffi+1) 


Moments  of  negative  order  exist  if  g is  large  enough.  If 


o > r,  2 ' s then 


T 


, -r.  ^ r(ct.4gj_r(a-r) 
^ r(a)  r(a-t^-r) 


(E.5) 


Y^"®  = r(a-fg)  r(g-r)  r(8-s) 
LU  ;u-a;  r(a+e-r-s) 


(E.7) 


and  the  harmonic  mean  option  will  be  available, 


1 + A: 


1 < a 


(E.8) 


Var(i) 


_ B(g+6-l) 


(g-1)^  (g-2) 


2 < g 


(E.9) 


and  if  g > 1,  B > 1 


(g-B-l) 
(g-1)  (B-1) 


(E.IO) 


Symmetric  Beta 

Set  g = B and  obtain 


’I 


I 


u = 1/2 


= 


4(2g+l) 


E(i)  = ^ 

g-1 


Var(|)  = -^L(2|=1I. 


(g-1)"  (g-2) 


^ ,1  1 . (2g-l) 

Cov(x,  - 2 

(a-1) 
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V -•«>  ■*  ■ m 


1 < g 


2 < g 


1 < g 


(E.ll) 

(E.12) 

(E.13) 

(E.14) 

(E.15) 


APPENDIX  F 


APL  PROGRAMS 


The  function  HARP  performs  the  iteration  (3.7)  to  estimate  the 
Poisson  parameter  from  the  harmonic  mean.  The  left  argument  X is  the 
initial  value  (usually  x)  and  the  right  argument  is  y of  (3.6).  The 
function  HMV  computes  the  variance  of  y using  (B.7)  and  (B.8).  The 
right  argument  is  the  set  of  parameter  values,  A. 


V L*-X  HARP  i'-,LL 

[1] 

[2]  L1:LL*-L 

[3] 

[4]  L 

[5]  ■*Ll»\i\L-LL)iO.OOQl 

V 


7 Y^HMV  u-.BiD  iMi'i 
Cl]  A’- 50 
[ 2 ] k-*-(  I.  J *1  ) “ . * I iV 

[ 3 ] £>-*-(  (.’•/•*- 1 tp  V)  ,/V)p(  *-iA)x(  \}])  x>  xH 

[4]  3*^1  H ,M)  a *-L 

[5] 

[6]  y-^YiL 

[ 7 ] r'*-  y~(  ( 1 - *-L  ) iB  ) * 2 

V 


The  function  MONT  produces  E(X  ) for  the  symmetric  beta  distri- 
bution using  (E.2)  with  a = 8.  The  left  argument  is  the  (integral)  order 
of  the  moment  and  the  right  argument  is  the  parameter.  The  function  VAR 
computes  the  variance  of  the  estimator,  a of  (3.16)  using  (3.17)  and 

(A. 2).  Again  the  argument  is  the  parameter 
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-% 


- « 


'Ik 


1 


V Z^R  MONT  A -.Nil  | 

Cl]  l-^o 

[2]  ZHN.N^pA)pl 

[3]  L1:J^J+1  I 

[4]  Zx-Zx(^(iy,;i?)p4+j.l  +4+J.1  I 

[5] -*’LlxiJ<^  j 

V i 


V V-*-VAR  A 

Cl]  ^>^(4  MONT  >l)-(4x(3  MONT  4 ) x(  l MONT  >l))-(3x(2  MONT  /I  ) *2  ) - 4x  ( 4+ 8x4  ) * -2 
C2]  7-H4X  Vx(l  + 2x.4  )*4 

V 


T\e  polygamma  functions  are  computed  using  PSI  and  JEX.  When  the 
(scalar)  left  argument,  N,  is  zero  the  psi  function  is  produced.  Integral 
values  of  N index  the  order  of  the  derivative  of  psi.  The  argument  of 
the  function  appears  on  the  right.  The  technique  comes  from  the  asymptotic 
expansions  in  Abramowitz. 


7 P*N  PSI  liCilViJIViX-.KKiYYiViZiTil 
Cl]  n N IS  THE  ORDER  OF  THE  DERIVATIVE  OF  THE  PSI  FUNCTION 

C2]  rt  y (>0)  IS  THE  ARGUMENT , SCALAR  OR  VECTOR 

C3]  C-^10 

C4]  IV*-\pY*,Y 

[5]  F*-Z-KHpY)pO 

[6]  KK*-fC-YZJIVHV*Y<C)/IVl 

C7]  -Ll»\  (pIV)=p(~V)/IV 

C8]  r-<-yc^iK] 

C9]  I-^O 

CIO]  L2;I-^Jtl 

Cii]  yy«^xxci]pTCi  ] 

Cl  2]  zCr]-^(  .'iV)x+/(  (yy-1  )+ixx[l])  *-!+// 

Cl3]  •*‘L2»\I<pJIV 
C14]  Z^V\Zl\pJIV] 

Cl  5]  K->rV\KK 
C16]  i,l;-*'51xi  JV>0 
Cl7]  F>--{9K*Y)-(2>^K+Y)*-1 
ClB]  -►2+12  6 

Cl9]  51:  P-(  ( .'//-I  )x(  y+X)  *-,V)+(  ;A^)xO.  5x(  Y+K)  *-N  + l 
C20]  P*-(.  Cl)*N  + l)»P+Z*N  JEX  Y+K 


I 


L 


. >■  -V. 


' » 


■’-f- 


V JEX  YiCiMiFiEiA 

Cl]  n USED  IN  THE  P0LyCA!4MA  FUNCTIONS 

[2  ] C-^-CH  5,(7t5),(5t7),(lli25),(  5*  455Ulx691  ),(6  91  i7x455)) 
[3]  ,i>)pC 

C 4 ] M*-l 

[5]  f-^(  y*2  )o  . x(  2+2x1  W-1  )x(  i+2x  iW-i  ) i(  /;r+2x  ) x(  ff+l  + 2x  iA/-i  ) 

[6]  E*-Cy.F 

[7]  7*6  )x(  y*-A/+2x,V)x{  ! ~l  + /7+2xW)  *(  ' 2*M) 

[8]  c7-<-4x1  + £’C;6]x1+£’C;53x1  +£"[  ;4]*l+£'[;3]xl+£’C  ;2]xl+£’[;l] 

V 


The  efficiency  of  the  usual  moment  estimator  (4.5)  of  gamma  distri- 
bution parameters  is  produced  by  EFF  using  (4.8)  and  (4.14).  The  efficiency 
of  (4.10)  is  computed  by  EFFH  using  (4.13)  and  (4.14),  and  the  relative 
efficiency  (4.15)  is  the  function  REFF.  Only  the  parameter  a is  needed 
and  it  is  entered  on  the  right  for  all  three.  It  must  be  >2  for  the 
latter  two  functions. 


V E*EFF  I 

[1]  £'<H(2x(y+i)x((y*(i  psj  y))-i))*-i 

V 


V E-^EFFH  y 

Cl]  £’-2*(y-2  )*(  y-1  )*2 

C2]  £'-H£x(yxi  PSI  y)-l 

C3]  £'-£•*-1 

7 


7 Z*-REFF  A 

Cl]  Z-^(4-1)x(4-1  )i(i4  + l )x/l-2 
7 
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#»  . 


We  turn  now  to  the  programs  that  support  the  negative  binomial 


distribution.  The  function  MM  computes  the  ordinary  method  of  moments 

estimators  p , r of  (4.25).  The  left  argument  is  x and  the  right 
2 

argument  is  s . 

V MM  S 

[1]  P^XiS 

[2]  E*(X*2)iS-X 

V 

The  determinant  of  the  asymptotic  covariance  matrix  of  p,  r is  produced 

(see  (4.29))  by  DM,  whose  left  and  right  arguments  are  r and  p. 

_ 2 

The  function  PSIB  computes  x (=  XB) , s (=  S) , y (=  Y) , the 
estimators  p,  r (using  MM),  and  the  maximum  likelihood  estimate  p,  r 
by  applying  Newton's  method  to  (4.18)  using  p,  f for  initialization. 

The  left  argument,  F,  is  the  vector  of  observed  frequencies  corresponding 
to  the  right  argument  J,  the  vector  of  variate  values. 


7 PSIB  J xPH’,PKD-,Z-,ZZ 

Cl]  XBH*/J*F)i*/F 

[2]  S*-iH  +/F)-1  )x(  +/FxtT*2)-(+/F)x;C5*2 

[3]  YH-*-/FiJ  -n)i-^/F 

[4]  XB  MM  S 

[5]  E,P 
[6j  Ll-.PP-^P 

[7]  2-^+/(0  PSI  EW)*Fi-¥/F 

[8]  ZZ-^+/(l  PSI  EW)*Fi-t/F 

[9]  PH*-Z-(.0  PSI  E)-{9E)~9E+XB 

[10]  PFO-t-ZZ-d  PSI  E)-iiE)-iEi-XB 

[11]  E->-E-PHiPHD 

[12]  F*-E*E+XB 

[13]  E,P 

[14]  -►L  lx  I ( I pp-p)20.  0001 

V 


The  function  INFI  computes  the  difficult  element  of  the  information  matrix, 
‘^22  (A. 20),  using  (D.26).  The  left  and  right  arguments  are  r and  p, 

resp.  The  vector  r must  be  whole  numbers  except  it  will  also  handle  the 
values  . 5,  1. 5,  ....  4. 5. 


V Z*-R  INFr 

Cl] 

[2] 

[3]  i4*(  (pA’)  ,pi?<-,P)pO 

[4] 

C 5 ] A 1 ! «/  ■♦*«/ 1 1 

[6]  /;+/?[«/])  >55 

[7]  -*L2*\  {aN)=QNN*(~V)/ K 

[a]  ZZ-^(  pA^l-*-F/A')pO 

[9]  i4C/Vl ; J / ( ZZo  . + x P[J]  -1  ) iNlo  . + i 

[10]  -►L3XI  (pAiV)=0 

[11]  A 2 : /I  [ A’-V  ; J-  ] -H(  : /VA' ) X ( .'  /f  [ J ] - 1 ) 5 I w [ j ] . 1 

[12]  A3:->‘AlxiJ<pP 

[13]  fl ( pP  ) , p .V  ) p /V* 2 

[14] 

7 


The  function  DETI  and  DI  both  compute  the  determinant  of  the  information 
matrix,  but  the  former  uses  INFI  in  (4.19)  and  can  handle  the  same  set 
of  r values  which  are  whole  numbers  plus  .5.  The  latter,  DI , uses  (D.25) 
and  all  values  of  r must  be  whole  numbers. 


T 


V D-*-R  DI  iZ 

[1] 

[2] 

[3]  Pl-^^(  (pP-^.P)  ,pP)pP*-2 

[4]  /1-h(  (pAf)  ,pP-K.i?)pO 

[5]  Z->-(p/|^)pO 

[6]  J*-Q 

[7]  L1:JW  + 1 

[8]  /lC;J]-x/(Z«.+iPC./]  )*P“.  + iPC.7] 

[9]  -*Ll*\J<pR 

[10]  (pif)  ,piV)piy+l 

[11]  0-ePl*D-^C+.  x/liS 

V 


The  determinant  (4.29)  of  the  asymptotic  covariance  matrix  (4.28) 
is  computed  by  the  function  DM 


7 P-^P  DM  P 

[1]  D-^2>«(  1+P*-.P)«.*  (P*2  )tl-P-H,P 

7 


and  the  efficiency  of  the  moment  estimator  is  NBEF 


7 E^R  REEF  P 

[1]  E*-i{R  DM  ?)x«P  DEri  P 
7 


The  function  HAR  implements  the  iteration  scheme  described  in 

* * 

(4.32)  and  produces  the  harmonic  mean  based  alternative  estimators  p , r 
from  the  system  (4.30).  The  left  argument  is  x and  the  right  argument 


. i 

I 


is  y . 


7 P-A  HAR  7;PP;C;C;CP;Z 


-CxP 


[1] 

P*Q->-Q.5 

[2] 

c->-y-i-Yxx 

[3] 

R<rX 

[4] 

Pi:  PF-hP 

[5] 

Z-P*R ) 

[6] 

£7P-*-C  + ZxZx- 

[7] 

P-^P-GiGP 

[3] 

R-X*PiQ->-l 

[9] 

R,P,G 

[10] 

-►P  1 X » ( 1 c ) : 
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The  Var(l/1+X)  is  computed  by  VHMl  using  (D.14)  and  (D.15)  for  r values 
that  are  whole  numbers,  and  using  (D.17)  and  (D.19)  for  r values  that 
are  half  way  between  whole  numbers. 


V Z*R  VHMl  PiAAiBBiBiIiPPiRRiSiHiHHiHHH;MiUliU2iJ;NiHNiQiSS 

[1]  fl  VAR  OF  i{l^X)  FOR  NEG  BIN(R,P)iR  MUST  BE  WHOLE  OR  WHOLE  PLUS 

[2]  RH~U2->-R  = Q.S)/ R*-(~Ul->-R  = l)/ R 

[3]  Z*-S'*-{  ( RR -*-p  R'*- ,R ) , PP'*-()P'*- ,P)p  0 

[4]  -►i[.3xi(pi?)=0 

[5]  I-HO 

[6]  L1:I^I+1 

[7]  B-Ril]-1 

[8]  -►(3+J26)xii3*L3 

[9]  B^3-l 

[10]  S[I;]*-aP 

[11]  -*■(  2+126  ) xiB=L5 

[12]  3[I;  ]-<-2*a2il  + P*0,  5 

[13]  BB*(PP.lB)pR[I]-li-iLB 

[14]  44-^(Po.  *-i?[I]-l+iLS)-l 

[15]  S[J;  ]*3[I;  ]+  + /4.UBB 

[16]  *LlxiZ<LRR<-pR 

[17]  L3:  P-^( -32) \P 

[16]  -►L4x;(p/?)  = 0 

[15]  R[U2/ I RR^RR++/ U23-G.5 

[20J  S^(~U2)\[1  ] S 

[21]  -o(  2+126  )xi  0=  + /i/2 

[ 22  ] S[  U2/ iRR;  ] -^+ ( 2 xa2  + 1 + P*0 . 5 ) +2  x(  ?*0.  5 ) -1 
[23  ] Z-^Jx(ff.Hi9Po  . *1?)  VJHH*-{  P-1  )o  . xi-p 

[24]  M*{(,HH*-{RR  ,PP)pP)-H)  iHHH 

[25]  Z*Z-M*2 

[26]  L 4 : /!/•»-(  2 X a 5 0 T Px  1 00 0 ) vaC"*”l -P 

[2  7]  A'[  ( V^x-,V<  50)/ i PPJ-^5  0 

[26]  SS  ■'•PPp  0 
[2  9]  .7-0 

[30]  i.  ^ : J "*~J  + 1 

[31]  ( Qo  . *1  A',v  ) i ( PP  ,iVA')p  ( \P/V-*-r  ECU  ] ) *2 

[32]  -»3  2 X \ u'  < PP 

[33]  -(  1+X26  )*( +/6T  ) *0 

[34]  R*-(~U1)\R 

[3  5]  P[6‘l  / i PP-PP  + + / O'l  ]*-l 

[36]  Z-^(~i/l)\[l  ] Z 

[3  7]  zrol/iPP;  ]-^Px(So-Px(  (ap)i,2)iQ)iQ 

V 


=1 


from  (4.36) . Two 


The  function  DMSI  computes  the  determinant  of  M 
auxiliary  functions  are  needed:  G22  provides  and  g^2  from  the 

second  row  of  A in  (4.33)  (compare  (2.17))  and  DMSI  is  needed  to  handle 
values  of  r = 1,  special  handling  being  required  becuase  of  (4.34). 


V DMSI  PiRRiPPiHH-,CiL-,GiG21;ZiU 

Cl]  AV( (pR)  ,pP)pO 

[2]  fl'«-(~i/-H/?  = l )//? 

C3]  Z^(.L*-R‘  .HP*2)^Q-*-l-P)^G->^R  G22  P 

[4]  ZH  Z+(«J(721  (/?i?-^p7?)  ,PP-^0P)pP)*2 

[5]  Z^ZJ(LxC-‘-(  ( (P-1  )»  . xC  )*2  )xP  VHMl  P)-«)G21*2 

[6]  .tC  (~i/)/ IPP++/P;  ]-«-Z 

[7]  -►(2  + l26)xi0  = +/i; 

[8]  ;rci//ipp++/p;]^p«si  p 

V 


V Z^DMSl  Piw 

Cl]  Z->-(  ( 0 . 5x(  «P  )*2  ) -p-ap)  *2 

C2]  ^ z^Z4(((pP)pi  mn  P)x((i-?)*3)).(vxp)*2 


V Z^R  G22  P;PP;PP;Z’;rp;P;PP;y;y 

Cl  ] H-*-(  (,PP-*-pP*-  ,P)  ,PP-*-pP-*-,P  )pPo  . 

C2]  nP-«>(PP.PP)pP 

[3] 

C4]  V^R  = X 

[5]  2’2’^(-/)/»PP 

C6]  Z-H^<5^{RR,PP)p9P 

[7]  Z[;i'PJ-^i’C;P2’]  + Z[;i’P]f(pZ[  ;P2’])p(~y)/p-i 

Cd]  C21-( (PP,PP)oP)xR4RR 

[5]  C2:-H(721-(l-P)jfi?(PP,Pp)pi-P 

V 


* A 

The  efficiency  of  the  harmonic  mean  based  estimator  p , r is  computed 
by  the  function  EFF 


7 E-R  EFF  P 

[1]  EHR  DMSI  P)i6lR  DETI  P 
V 


Both  EFf  and  NBEF  used  DETI  and,  for  that  reason,  are  restricted  to  only 

■k  * 

a few  fractional  values  of  r.  The  relative  efficiency  of  p , r with 
respect  to  p,  r is  not  so  restricted.  It  is  computed  by  RELEF  and  accepts 
any  r > 0. 

7 E->-R  RELEF  P 

Cl]  E*-{{R  DM  P)»{R  DMEI  P)) 

7 


Methods  for  computing  the  efficiency  of  beta  distribution  estimators 
are  supported  by  the  following  programs:  The  determinant  of  the  information 

matrix  (4.39)  is  computed  by  the  functions  DINF  and  DDINF.  The  former  takes 
a single  (vector)  argument  and  produces  a symmetric,  sq’.are  matrix  of 
values  |a|  for  all  pairs  of  components  of  the  arguments.  If  the  arguments 
a and  3 must  be  ent  red  separately,  then  the  two  argument  DDINF  can  be 
used. 


7 u->-DINF  AiiViBiC 
Cl]  L-<-  ,L*-A  ° . +A 
C2]  5-^( /V  ,;/■<- p/1  ) 0 1 PST  L 

C3]  I.-(Co  . xC  )- (Co  . FSI  A)*B 

7 


V L -A  DDINF  B;M-,N  iCA  iCBiD 


Cl] 

L*-  ,L  *-A  o . + 3 

C2] 

D-^(.  {M*-oA  ) ,N*-pB  )p  1 

PS  I L 

C3] 

L-(C/lo  ,^CB)-{  (.CA*-1 

PS  I A ) » . + C5-^1  PS  I B)>^D 

t 7 
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^ ^ 


j 


L 
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The  function  DETMM  accepts  two  by  two  matrices  A (on  the  left) 

2 

and  C (on  the  right)  and  computed  the  determinant  |m|  = |c|/|a| 
of  (2.7).  The  function  ADET  uses  it  to  produce  an  array  of  such 
determinants.  The  arguments  H and  C are  four-dimensional  and  may  be 


thought  of  as  an  M by  N array  of  2 by  2 matrices.  The  left  set,  H, 
are  the  coefficient  arrays  A of  (2.2)  and  the  right  set,  C,  are  the 


covariances  (2.6) 


V M->-A  DET14M  C 
Cl]  + 

C 2 ] <V-^(  /■^C  1 ; 1 ]x.V[2  ; 2 ] ) -W[  1 ; 2 ] *2 

7 


V MM*-K  ADET  CiN\I-,J\M 

[1]  2 + pfl)  )p0 

[2] 

[3] 

C 4 ] C 

[5] 

[6]  iliJl  DETMM  C’[;;I;J] 

[ 7 ] -*L2*xJ  <N 

[ 8 ] 1 X 1 1 <i>/ 

V 


The  computation  of  the  efficiency  of  the  ordinary  moment  estimator 


(4.42)  requires  the  coefficient  array  (4.43)  and  the  convariance  array 
(4.45).  The  former  is  computed  by  the  function  COEFM  and  the  latter  by 
COVM.  Each  takes  a single  vector  argument  and  computes  the  required 
values  for  all  pairs  of  components  in  the  argument.  The  function  COVM 
requires  the  beta  distribution  moments  (E.2)  and  these  are  computed  by 
MONT. 
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V H-^COEFM  A-,N‘,B 
Cl]  //*(  2, 2,/V.iV-p/4  )pl 
[2]  //Cl;l;;]-^(iV.JV)p-4 

C3]  //Cl;  2;  ; ]-HS)(//,A?)p/t 

C4]  S*(/l<>.x<4)x(j4®.+.4)i(i4'».+.4+l) 

C5]  //C2;l;  ; + -i4  )x(//,.V)p^ 

C6]  ffC2;  2;  ; ]-t-5+(  ( --4  ) <»  . +/1  ) x«j  ( A' , // ) P>1 

C7]  H^HH2,2,N  ,N)pAo  .-^A 

1 


7 Z->-R  MONT  AiNiI 
Cl]  I^Q 

C2]  Z-^-C  A )pi 

C3]  Ll:I-<-J+l 

C4]  2-<-Zx(«)(  A ,A)p4+I-l  )ii4®.  +4+1-1 

C5j  I <R 

7 


All  these  are  utilized  by  EFBM  which  computes  the  efficiency  (2.22). 

The  output  is  a symmetric  matrix.  The  argument  must  have  positive  components. 


7 E^EFBM  AiLiC-.HiM 

Cl]  L*-DINF  A 

C2]  C^CCVM  A. 

C3]  H^COEFM  A 

C4]  ADET  C 

C5]  E-*-iM'*-L 

7 


« 

=» 


I 


•« 


The  efficiency  of  (4.48)  is  handled  in  similar  fashion.  (Also 
with  single  arguments.)  An  array  of  2 by  2 matrices  (4.49)  is  produced 
by  the  function  COEF  and  the  matching  matrices  (4.50)  by  COV.  These  are 
used  by  EFBH  to  compute  a syrometric  matrix  of  efficiencies.  The  argument 
must  have  all  components  >2.  1 

i 

jl 
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V H-^COEF  A‘,N 

Cl]  HH2,2.N  ,N^pA)p-l 

[2]  //C2;2;  ;]<^/l‘>.M-l 

[3]  ffCl;l;  ;]^«?ffC2;2;  ;] 

V 


V Ct-COV  A;  A’ 

Cl]  2. 2,A,A-Hp4  )pl 

C2]  CC  1;  1 ; ; ]-*-(A  9 . +A -1  ) ® . M -2 

C3]  CC1;2;  ;]«^CC2;1;  ;]<h-46.+4-i 

C4]  CC2;2;  ;]-<-{A«».+4-l  )x(A<».  M-2  ) 

V 


V E-EFBH  A iCiHiL 

Cl]  L->-DIEF  A 

C2]  H^COEF  A 

C3]  C*-COV  A 

C4]  M-ff  ADET  C 

C5]  E*iM»L 

V 


The  estimator  (4.52)  is  managed  in  like  fashion,  only  this  time 
the  arguments  a (left)  and  6 (right)  must  be  entered  separately  with 
a > 2,  6 > 0.  The  coefficients  (4.53)  are  computed  by  the  function  COEFH 
and  the  covariances  (4.54)  by  COVH.  These  are  drawn  on  by  EFBMH  to  produce 
the  efficienices 


V FI->-A  COEFH 

Cl]  H*-{2,2.{M*-pA)  ,N*-pB)p-L 

C2]  A?Cl  ;1;  ;]-^5(Ao.+5)t(W,,V)p-5 

C3]  i?Cl;2;  ;]-H((5)(yv,,V)D4)i/lo.+B 

C4]  nC2;l;  ;]-<-J(/l-l  )o  . iS 


. -hiarr  t u — ■ ■ *^i  n» , r . v 


V C--A  COVH  BiNiM 

Cl]  C*{2,2,(M^pA)  ,N->-pB)pO 

[2]  CCl;l  ;;]<*-(/l«.*fl)i/lo.+s  + i 

[3]  <7Cl;2;  ;J-K7C2;1;  il-(M,N)p-B 

[4]  (7C2;2;  ;]-^(i4».+S-l)t(il-2)«.»fl 

V 


V E*-A  EFBMH  BiC-,HiM\L 

[1]  L-^A  DDINF  B 

[2]  H*-A  COEFH  B 

[3]  C+-A  COVH  B 

[4] '  M*-H  ADET  C 

[5]  E*-iM*L 
7 
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